# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from matplotlib.pyplot import acorr
from scipy.stats import t, shapiro, kstest, uniform, anderson, norm
from algorithm_A import perform_algorithm_A
from sklearn.preprocessing import StandardScaler
from perform_anderson_test import perform_anderson
from autocorrelation_test import perform_autocorrelation_test
from statsmodels.tsa.stattools import adfuller
from linear_regression_test import perform_regression_test


def calculate_risk_1(x):
    risk_1 = 0
    x, x_star, s_star, higher_bound, lower_bound = perform_algorithm_A(x)
    flag = 0 
    for i in range(len(x)):
        if ~flag:
            n = 0
        if x[i] == higher_bound or x[i] == lower_bound:
            flag = 1
            n += 1
            risk_1 += n
        else:
            flag = 0
        
    risk_1 /= len(x)
    
    # 判定数据是否像传感器 8 和 35 一样的平滑
    flag_robust = False
    if abs(lower_bound-higher_bound) <= lower_bound*0.01:
        #数据类似传感器 8 一样的，全部等于一个数的序列，所以置零。
        risk_1 = 0
        flag_robust = True
    
    return x, risk_1, flag_robust

# ks test 均匀分布检验需要注意的点：https://stackoverflow.com/questions/25208421/how-to-test-for-uniformity
def calculate_risk_2(x, flag_robust):
    '''
    根据 AD 检验和 KS 检验的 p-value 计算非独立性和偶发性风险（见博客）
    
    x 是经过算法 A 过滤后的值
    '''

    x_mean = np.mean(x)
    x_std = np.std(x)
    
    x_q3 = np.quantile(x, q=.75)
    x_q1 = np.quantile(x, q=.25)
    
    # AD 检验分析白噪声误差（正态分布）
    # p_ad = perform_anderson(x_norm)
    
    # 检验分析白噪声误差（正态分布)，若因为 0 除错误无法计算 p-value，则设为 0
    # 即让风险最大（下同）
    _, p_ad = kstest(x, norm(loc=x_mean, scale=x_std).cdf)
    if p_ad != p_ad:
        p_ad = 0
    
    # KS 检验分析平稳性（均匀分布）
    x_between = x[(x<=x_q3) & (x>=x_q1)]
    delta = abs(x_q3-x_q1)
    _ , p_ks = kstest(x_between, uniform(loc=x_q1, scale=delta).cdf)

    if p_ks != p_ks:
        p_ks = 0
        
    # 是否通过 ADF 检验
    flag = False
    # ADF 检验
    if ~flag_robust:
        res = adfuller(x)
        p_adf = res[1]
        if p_adf <= 0.05:
            flag = True
    else:
        p_adf = 0
    
    
    # 若 时间序列的标准化后，方差比较小，则将风险置 0 
    # flag = flag or perform_regression_test(x)
    
    risk_2 = 1-max(p_ad, p_ks)
    risk_2 = min(risk_2, p_adf)
    
    return risk_2, flag

def calculate_risk_3(x, flag, flag_robust):
    '''
    计算持续性、联动性风险
    '''
    series = pd.Series(x)
    risk_3 = 0
    if flag_robust or flag:
        return 0
    for lag in range(1, 9):
        # 自相关系数
        r = series.autocorr(lag=lag)
        if r != r:
            # 若 r 因为 0 除计算不出来
            r = 0
        
        n = series.shape[0] - lag
        try:
            r = perform_autocorrelation_test(r, n)
        except:
            pass
        # 加权和
        risk_3 += r*lag/sum([1,2,3,4,5,6,7,8])
    
    return risk_3


def calculate_risk(data):
    '''
    计算风险
    '''
    # 删除第一列（时间编号）
    data.drop(data.columns[0], axis=1, inplace=True)
    cols = data.columns
    
    risk_table = {}
    for col in cols:
        x = data[col]
        x = x.values
        x_after, risk_1, flag_robust = calculate_risk_1(x)
        risk_2, flag = calculate_risk_2(x_after, flag_robust)
        risk_3 = calculate_risk_3(x, flag, flag_robust)
        
        print('第一类风险为： ', risk_1, '\t第二类风险为： ', risk_2, '\t第三类风险为: ', risk_3)
        
        risk_table[f'传感器{col}'] = [risk_1, risk_2, risk_3, risk_1+risk_2+risk_3]

    table = pd.DataFrame.from_dict(risk_table, orient='index')
    table.columns = ['风险1', '风险2', '风险3', '总风险']
    table.to_excel(r'../附件/中间数据/风险评估表.xlsx')

if __name__ == '__main__':
    path = r'../附件/附件1(Appendix 1)2021-51MCM-Problem C.xlsx'
    data = pd.read_excel(path, index_col=1, header=0, parse_dates=True)

    calculate_risk(data)






